!========================================================================================
  program Test_Radiation 
!========================================================================================

  use HDF5
  use radiation_m
  use dtm_m
  use param_m
  use topology_m

  implicit none

  character(255), parameter  :: Test_radid   ='$Id: Test_Radiation.f90 34 2013-07-23 21:43:05Z vlecous1 $'
  character(255), parameter  :: Test_radrev  ='$Revision: 34 $'
  character(255), parameter  :: Test_raddate ='$Date: 2013-07-23 21:43:05 +0000 (Tue, 23 Jul 2013) $'

  real, dimension(:,:,:),   allocatable :: temp, pressure, div_qr
  real, dimension(:,:,:,:), allocatable :: yspecies

!-----------For code performance assesment only------------------------------------------
  real :: start       ! For CPU_TIME purposes
  real :: finish      ! For CPU_TIME purposes
  real :: init_DTM    ! For CPU_TIME purposes
  real :: Call_DTM    ! For CPU_TIME purposes
  real :: total       ! total = finish - start
!----------------------------------------------------------------------------------------

  integer :: io ! Unit of output file
  integer :: i_config ! Unit of the configuration input file (READ ONLY)

  integer :: npx, npy, npz ! Number of processors on x, y, and z
!----------------------------------------------------------------------------------------
  io = 6

  call CPU_TIME(start)

!----------------------------------------------------------------------------------------
! Open Configuration.in.
! Get number of processors in x, y, and z directions: npx, npy, npz
! Get total number of points in x, y, and z

  i_config = 2

  open(i_config,file='Configuration.in',status='old',form='formatted',action='read')
  read(i_config,*) 
  read(i_config,*) npx
  read(i_config,*) npy
  read(i_config,*) npz
  read(i_config,*) 
  read(i_config,*) nx_g
  read(i_config,*) ny_g
  read(i_config,*) nz_g
  read(i_config,*) 
  read(i_config,*) xmin
  read(i_config,*) xmax
  read(i_config,*) ymin
  read(i_config,*) ymax
  read(i_config,*) zmin
  read(i_config,*) zmax
  read(i_config,*) 
  read(i_config,*) periodic_x
  read(i_config,*) periodic_y
  read(i_config,*) periodic_z
  close(i_config)
!----------------------------------------------------------------------------------------

  open(io,file='Radiation3d.dat')

  numdim = 3

  call initialize_topology(io,npx,npy,npz)
  call Init_MPI_Rad
  call generate_grid
  call initialize_radiation(io)

  call CPU_TIME(init_DTM)
!----------------------------------------------------------------------------------------
! Create Temp, pressure, yspecies arrays

! Allocate memory for temp, pressure, div_qr, and yspecies

  if (allocated(temp))     deallocate(temp)
  if (allocated(pressure)) deallocate(pressure)
  if (allocated(div_qr))   deallocate(div_qr)
  if (allocated(yspecies)) deallocate(yspecies)

  allocate(temp(nx,ny,nz))
  allocate(pressure(nx,ny,nz))
  allocate(div_qr(nx,ny,nz))
  allocate(yspecies(nx,ny,nz,nsc+1))

  if (myid == 0) then
    write(io,*) 'Memory allocated for Temperature, Pressure, div_qr & yspecies.'
  endif

  call init_field(temp,pressure,div_qr,yspecies,io)

!----------------------------------------------------------------------------------------

  call DTM (temp,pressure,yspecies,div_qr,io)

  call CPU_TIME(Call_DTM)

  call write_savefile(temp, pressure,div_qr,io)

  call Print_irrad(myid)

  call Print_Fluxes(myid,io)

!  call Test_rays()

!  call Print_all_rays()

  call Print_Cell_Source(myid)

!  call Print_Diag_North

!  call Print_Div_qr(div_qr)
  
  call allocate_radiation_arrays(-1)

  call CPU_TIME(finish)

  total = finish-start

  if (myid==0) then
  print *, ' '
  print *, '-------------------------------------------------------------------'
  print *, 'Time| Total:                 ', total, ' (s).'
  print *, 'Time| Initialisation of DTM: ', init_DTM-start, ' (s).'
  print *, 'Time| Run of DTM:            ', Call_DTM-init_DTM , ' (s).'
  print *, '-------------------------------------------------------------------'
  print *, ' '
  print *, 'Program stopped normally. Goodbye.'
  print *, '-------------------------------------------------------------------'

  close(io)

!---FORMAT STATEMENTS--------------------------------------------------------------------
10 format('points number: ',I4,', ray number: ',I4)
11 format('number of cells crossed: ', I4)
12 format('Origin of the ray: x= ', 1f16.10, ' y= ', 1f16.10,' z= ', 1f16.10)
13 format('There is ',I4,' problems.')

100 format('index in x : ',I4,' in y: ',I4,'in z: ',I4,' length: ',1f8.4) 
 
  endif
!----------------------------------------------------------------------------------------
  end program Test_Radiation

!========================================================================================
  subroutine generate_grid
!========================================================================================

  use param_m
  use radiation_m, only : myid_cart, cart_dim 
  implicit none

  real :: x_min, y_min, z_min 
  real :: x_max, y_max, z_max
  real :: delta_x, delta_y, delta_z
  real :: xid, yid, zid  

  integer i

!----------------------------------------------------------------------------------------
! Define structure of the grid

  nsc = 1 

  l_ref = 1.0e-2 

  xid = myid_cart(1)
  yid = myid_cart(2)
  zid = myid_cart(3)

  x_min = xmin
  y_min = ymin
  z_min = zmin

  x_max = xmax
  y_max = ymax
  z_max = zmax

  delta_x = (x_max-x_min)/float(nx_g-1)
  delta_y = (y_max-y_min)/float(ny_g-1)
  delta_z = (z_max-z_min)/float(nz_g-1)

  nx = nx_g/cart_dim(1)
  ny = ny_g/cart_dim(2)
  nz = nz_g/cart_dim(3)

  allocate(x(nx))
  allocate(y(ny))
  allocate(z(nz))

  do i = 1, nx
    x(i) = x_min + float(i-1)*delta_x + nx*delta_x*xid
  enddo

  do i = 1, ny
    y(i) = y_min + float(i-1)*delta_y + ny*delta_y*yid
  enddo

  do i = 1, nz
    z(i) = z_min + float(i-1)*delta_z + nz*delta_z*zid
  enddo

! Non dimensionlize

  x = x/l_ref
  y = y/l_ref
  z = z/l_ref

!----------------------------------------------------------------------------------------
  return
  end subroutine generate_grid


!========================================================================================
  subroutine init_field(temp,pressure,div_qr,yspecies,io)
!========================================================================================
! This subroutine initialize the Temperature, Pressure, yspecies, and Divergence of Qr
! Dimensions:
!------------ 
!  Temp(nx,ny,nz)
!  Pressure(nx,ny,nz)
!  yspecies(nx,ny,nz,nsc+1)
!  Div_qr(nx,ny,nz)
!----------------------------------------------------------------------------------------
  use topology_m,  only : myid
  use radiation_m, only : myid_cart, cart_dim 
  use radiation_m, only : Temp_South, Temp_North, Temp_Right, Temp_Left, Temp_Bottom, Temp_Top
  use param_m,     only : nx, ny, nz, nsc, x, y, z, l_ref

  implicit none

! Variables IN and OUT
  real,intent(inout), dimension(nx,ny,nz) :: temp
  real,intent(inout), dimension(nx,ny,nz) :: pressure
  real,intent(inout), dimension(nx,ny,nz) :: div_qr

  real,intent(inout), dimension(nx,ny,nz,nsc+1) :: yspecies

  integer, intent(inout) :: io

! Local Variables
  real :: xid, yid, zid ! Cartesian coordinates of the local processor
  real :: x_s, y_s, z_s, radius_sphere, radius ! Variable for hot sphere test
  
  integer :: j_x, j_y, j_z


! Define configuration of the hot sphere
  x_s = 0.0d+0
  y_s = 0.0d+0
  z_s = 0.0d+0

  radius_sphere = 0.3d+0 ! meters

  xid = myid_cart(1)
  yid = myid_cart(2)
  zid = myid_cart(3)

  temp     = 1.27315d+3    ! Units in Kelvin

! Compute distance from center of the sphere
! If the point is inside the sphere, assign it temperature to be 1200 K
!----------------------------------------------------------------------------------------

!   do j_z = 1, nz
!     do j_y = 1, ny
!       do j_x = 1, nx
! 
!         radius = dsqrt( (x(j_x)*l_ref-x_s)**2 + (y(j_y)*l_ref-y_s)**2 + (z(j_z)*l_ref-z_s)**2  )
!         if (radius <= radius_sphere) temp(j_x, j_y, j_z) = 2.0d+3
! 
!       enddo
!     enddo
!   enddo


  pressure = 101325.0d+0  ! Units in Pascal
 
! Define Faces temperatures
!----------------------------

  Temp_South  = 1.0d-5
  Temp_North  = 1.0d-5
  Temp_Right  = 1.0d-5
  Temp_Left   = 1.0d-5
  Temp_Bottom = 1.0d-5
  Temp_Top    = 1.0d-5

  div_qr = 0.0d+0

  yspecies(:,:,:,1) = 0.233d+0
  yspecies(:,:,:,2) = 1.0d+0-0.233d+0


  if (myid == 0) then
    write(io,*) 'Fields initialized.'
  endif

!----------------------------------------------------------------------------------------
  return
  end subroutine init_field